Life History Strategies in Lizards: Insights into the Slow-Fast Continuum

Set Up

Initilisation

Load Libraries

library(tidyverse)
library(corrplot) # For correlation plots
library(FactoMineR) # For PCA
library(factoextra) # For PCA plots
library(vegan) # For CCA
library(ggplot2)
library(rsample) # For sampling of data (split, analysis,..)
library(rpart)
library(rpart.plot)
library(randomForest)
library(tibble)
library(tidyr)
library(gridExtra)
library(factoextra)
library(caret)
library(cluster)

Knitting Configuration

# Set to render chunks during the knitting process
knitr::opts_chunk$set(echo = TRUE)

Set Up the Environment

# Clear Environment
rm(list = ls())

# Set Working Directory
setwd("/Users/eamon/Desktop/MathToolsProject/R_Project/MathematicalTools/")

# Load Data
data <- read.csv("lizard.csv")

View the Data

str(data)

'data.frame':   740 obs. of  24 variables:
 $ Species                        : chr  "Agama Agama" "Agama Agama" "Agama hispida" "Caimanops amphiboluroides" ...
 $ Genus                          : chr  "Agama" "Agama" "Agama" "Caimanops" ...
 $ Family                         : chr  "Agamidae" "Agamidae" "Agamidae" "Agamidae" ...
 $ Population                     : chr  "Suacocco, Liberia 1952-1953" "Ibadan, NIgeria" "Kalahari Desert, Africa, 1969-1970" "Australia" ...
 $ Longitude                      : num  -9.05 3.92 21.25 119.1 119.62 ...
 $ Latitude                       : num  7.08 7.4 -27.03 -28.5 -27.08 ...
 $ Average.Female.adult.weight..g.: num  NA NA 31.9 NA 10 ...
 $ SD.Female.adult.weight..g.     : chr  NA NA "6.16" NA ...
 $ Sample.Size.Female.adult.weight: int  NA NA 28 NA 17 75 NA 41 741 120 ...
 $ Mean.F.SVL.adults..mm.         : num  85 97 94.8 94 62.9 ...
 $ SD.F.SVL.adults..mm.           : num  NA NA 6.91 7.07 4.43 ...
 $ Sample.size.Mean.F.SVL.adults  : int  63 NA 51 2 18 75 NA 41 756 120 ...
 $ F.SVL.at.maturity..mm.         : num  75 NA 76 89 57 44 48 49 47 51 ...
 $ Offspring.SVL..mm.             : num  32 NA 29 38 22 23 24 25 22 26 ...
 $ Mean.Clutch.Size               : num  7.5 5.5 13.4 12 5.36 ...
 $ Sample.size.Clutch.Size.       : int  NA NA 45 2 14 28 NA 13 40 39 ...
 $ Mode.of.reproduction           : chr  "Oviparous" "Oviparous" "Oviparous" "Oviparous" ...
 $ Clutches.per.year              : num  2 NA NA NA NA NA 2 NA NA NA ...
 $ Clutch.Frequency               : num  2 2 1 1 1 1 2 1 2 2 ...
 $ RCM                            : num  NA NA 0.305 NA 0.166 ...
 $ Foraging.Mode                  : chr  "Sit and wait" "Sit and wait" "Sit and wait" "Sit and wait" ...
 $ Distribution                   : chr  "Tropical" "Tropical" "Temperate" "Temperate" ...
 $ Prefered.Habitat.Type          : chr  "Saxicolous" "Saxicolous" "Semi-arboreal" "Semi-arboreal" ...
 $ Source                         : chr  "(Harris 1964)" "(Daniel 1960)" "(Pianka 1986)" "(Pianka 2013d)" ...

Cleaning

Change Variable Types

# This variables was input as a character, but it should be numeric
data$SD.Female.adult.weight..g. <- as.numeric(data$SD.Female.adult.weight..g.)

Warning: NAs introduced by coercion

Set Blanks in Categorical Variables as NA

# Convert the blank values of categorical variables to NA
data<-data %>%
  mutate(across(where(is.character), ~ na_if(.x, ""))) %>%
  mutate(across(where(is.character), as.factor))

Check for and Remove Duplicates

# Save duplicates as an object
duplicates <- data[duplicated(data), ]

# Remove duplicates
data <- data[!duplicated(data), ]

Visualise the Proportion of Each Variable Which is NA

By doing this we found that most of the numeric variables had a high proportion of NA’s, with some such as RCM having over 50%. Therefore exclusion of all rows containing NA’s was likely to drastically reduce the size dataset and cause data loss. The number of rows in the original dataset was 734 and the number after exclusion of all NAs was 146. Therefore the NA’s were removed separately for each analysis using only the variables required for that analysis, which helped to reduce data loss.

# Calculate percentage of NA values for each column
na_percentage <- colSums(is.na(data)) / nrow(data) * 100

# Convert NA percentages to a data frame
na_df <- data.frame(
  Column = names(na_percentage),
  NA_Percentage = na_percentage
)

# Create bar plot
ggplot(na_df, aes(x = Column, y = NA_Percentage)) +
  geom_bar(stat = "identity", fill = "skyblue") +
  labs(
    title = "Percentage of NA Values per Column",
    x = "Columns",
    y = "Percentage (%)"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1) # Rotate x-axis labels
  )

# Print the number of rows in the full dataset
print(paste("Number of rows in full dataset =", nrow(data)))

[1] "Number of rows in full dataset = 734"

# Print the number of rows after NA removal
print(paste("Number of rows after NA removal on entire dataset =", nrow(na.omit(data))))

[1] "Number of rows after NA removal on entire dataset = 146"

PCA

Choose the Variables to be Used in the PCA

This saves the variables to be used in the PCA as an object. The variables used are easily altered by commenting them in or out. Only numerical variables are considered as candidates as PCA can only use numerical variables. Where two variables describe the same trait and/or are correlated, one of them was removed. This was done in order to avoid bias in the PCA by over emphasising particular traits.

“F.SVL.at.maturity..mm” describes the same trait as “Mean.F.SVL.adults..mm.”. A correlation of these two variables revealed a 97% correlation. Therefore, “F.SVL.at.maturity..mm.” was randomly chosen to be removed.

Clutch Frequency and Clutches Per Year were also correlated and describe the same trait and therefore Clutch Frequency was removed.

Sample sizes and standard deviations were removed from the analysis as they were not thought to be relevant.

variables_to_use_PCA <- c(
  ### Location
  #"Longitude",
  #"Latitude",
  
  ### Means and Numerical
  #"Average.Female.adult.weight..g.",
  "Mean.F.SVL.adults..mm.",
  #"F.SVL.at.maturity..mm.",
  "Offspring.SVL..mm.",
  "Mean.Clutch.Size",
  "Clutches.per.year",
  #"Clutch.Frequency",
  "RCM"
  
  ### Standard Deviations and Sample Sizes
  #"SD.Female.adult.weight..g.",
  #"SD.F.SVL.adults..mm.",
  #"Sample.Size.Female.adult.weight",
  #"Sample.size.Mean.F.SVL.adults",
)

Investigation of the Similar Variables

To investigate the correlations between variables.

Create New Dataset retaining only the PCA variables

Here a new subset of the data is created for the PCA. The PCA data was cleaned of all rows containing NA’s and standardised such that each variable has mean 0 and variance 1 (which allows for comparison of variables on different scales in the PCA). A second data set was created containing the rows retained after NA omition of the PCA data, but containing all the original variables. This is later used to color the PCA visualisations with categorical variables.

# Create the new PCA dataset
data_PCA <- data %>%
  select(any_of(variables_to_use_PCA)) %>% # filter variables
  na.omit() %>% # omit NA's
  mutate_all(.funs = scale) %>% # standardise the data
  rownames_to_column(var = "row_id") # set rownames as a column

# Create ancillary dataset with all original variables but filtered to keep only rows present in the PCA dataset
# This will be used for coloration of the PCA graphs with categorical variables

data_PCA_all_vars <- data %>%
  rownames_to_column(var = "row_id") %>% # set rownames as a column
  filter(row_id %in% data_PCA$row_id) %>% # filter the data to match rows of data_PCA
  column_to_rownames(var = "row_id") # Convert column back to rownames

# Convert the column back to rownames for the PCA data
data_PCA <- data_PCA %>%
  column_to_rownames(var = "row_id")

Visualise with a Correlation

This shows the correlation of each variable pairs in filtered PCA data. None of the variables are highly correlated with each other.

### CorrPlot
data_PCA %>% 
  cor(use="pairwise.complete.obs") %>% # Calculate the empirical correlation matrix
  corrplot() # Then graph this matrix

Run the PCA

Here the PCA is run and the results are presented in various forms including the output tables, a scree plot showing the percentage of explained variance by each principle component and the correlation of each variable with each principle component.

# Run the PCA
result_pca <- PCA(data_PCA, 
               scale.unit = TRUE, # Option to center and scale data (useless here)
               ncp = 18, # Number of components to keep (here, all)
               graph = FALSE)

# Results of the PCA
result_pca$eig

       eigenvalue percentage of variance cumulative percentage of variance
comp 1  1.9331037              38.662073                          38.66207
comp 2  1.6070287              32.140575                          70.80265
comp 3  0.7868455              15.736911                          86.53956
comp 4  0.4494887               8.989773                          95.52933
comp 5  0.2235334               4.470668                         100.00000

result_pca$var

$coord
                           Dim.1      Dim.2       Dim.3       Dim.4       Dim.5
Mean.F.SVL.adults..mm. 0.6976592  0.6159273  0.04570107 -0.18884879 -0.31008511
Offspring.SVL..mm.     0.3783690  0.8203961  0.18805831  0.29955692  0.24225387
Mean.Clutch.Size       0.8383756 -0.2943332 -0.13505408 -0.36570314  0.24190071
Clutches.per.year      0.3758802 -0.5275783  0.75400593  0.09854526 -0.04625082
RCM                    0.6775343 -0.4354788 -0.40326973  0.42501824 -0.08965820

$cor
                           Dim.1      Dim.2       Dim.3       Dim.4       Dim.5
Mean.F.SVL.adults..mm. 0.6976592  0.6159273  0.04570107 -0.18884879 -0.31008511
Offspring.SVL..mm.     0.3783690  0.8203961  0.18805831  0.29955692  0.24225387
Mean.Clutch.Size       0.8383756 -0.2943332 -0.13505408 -0.36570314  0.24190071
Clutches.per.year      0.3758802 -0.5275783  0.75400593  0.09854526 -0.04625082
RCM                    0.6775343 -0.4354788 -0.40326973  0.42501824 -0.08965820

$cos2
                           Dim.1      Dim.2       Dim.3       Dim.4       Dim.5
Mean.F.SVL.adults..mm. 0.4867283 0.37936643 0.002088587 0.035663864 0.096152777
Offspring.SVL..mm.     0.1431631 0.67304969 0.035365929 0.089734347 0.058686935
Mean.Clutch.Size       0.7028736 0.08663202 0.018239606 0.133738786 0.058515951
Clutches.per.year      0.1412859 0.27833883 0.568524941 0.009711168 0.002139139
RCM                    0.4590527 0.18964177 0.162626476 0.180640503 0.008038593

$contrib
                           Dim.1    Dim.2      Dim.3     Dim.4     Dim.5
Mean.F.SVL.adults..mm. 25.178595 23.60670  0.2654381  7.934319 43.014950
Offspring.SVL..mm.      7.405868 41.88162  4.4946469 19.963651 26.254214
Mean.Clutch.Size       36.359852  5.39082  2.3180669 29.753539 26.177722
Clutches.per.year       7.308761 17.32009 72.2536906  2.160492  0.956966
RCM                    23.746924 11.80077 20.6681576 40.187999  3.596149

# Plot the Scree Plot to Show Variance Explained by Each Dimension
fviz_eig(result_pca, choice = "variance")

# Plot the Correlation of Each Variable and the Dimension
corr_plot_pca <-corrplot(result_pca$var$cor)

corr_plot_pca

$corr
                           Dim.1      Dim.2       Dim.3       Dim.4       Dim.5
Mean.F.SVL.adults..mm. 0.6976592  0.6159273  0.04570107 -0.18884879 -0.31008511
Offspring.SVL..mm.     0.3783690  0.8203961  0.18805831  0.29955692  0.24225387
Mean.Clutch.Size       0.8383756 -0.2943332 -0.13505408 -0.36570314  0.24190071
Clutches.per.year      0.3758802 -0.5275783  0.75400593  0.09854526 -0.04625082
RCM                    0.6775343 -0.4354788 -0.40326973  0.42501824 -0.08965820

$corrPos

$arg
$arg$type
[1] "full"

Plot the Individuals with the Variable Arrows and the Categorical Coloring

By altering the object Color_label we were able to check the effect of coloring with each of the four main categorical variables (Distribution, Habitat type, Foraging mode and Reproductive type). It was found that distribution produced the clearest clustering with the least noise. Therefore distribution was taken forward as a potential categorical variable which can account for much of the variation. This code could of been shortened using a loop, however it affected the scaling of the axes and the font size and we were unable to resolve this issue.

# Define Graphing Function
plot_pca_with_color <- function(color_label) {
  p1 <- fviz_pca_biplot(result_pca, axes = c(1, 2), col.ind = data_PCA_all_vars[[color_label]], col.var = "black")
  p2 <- fviz_pca_biplot(result_pca, axes = c(1, 3), col.ind = data_PCA_all_vars[[color_label]], col.var = "black")
  p3 <- fviz_pca_biplot(result_pca, axes = c(2, 3), col.ind = data_PCA_all_vars[[color_label]], col.var = "black")
  p4 <- fviz_pca_biplot(result_pca, axes = c(3, 4), col.ind = data_PCA_all_vars[[color_label]], col.var = "black")
  grid.arrange(p1, p2, p3, p4, nrow = 2, ncol = 2)
}

# Use function with Distribution
dist_plot <- plot_pca_with_color("Distribution")

# Save Plot
ggsave("PCA_plot_distribution.png", plot = dist_plot, width = 10, height = 8.4, units = "in", dpi = 500)
# Print plot
dist_plot

TableGrob (2 x 2) "arrange": 4 grobs

# Plot for other categorical variables
dist_plot <- plot_pca_with_color("Prefered.Habitat.Type")

dist_plot <- plot_pca_with_color("Mode.of.reproduction")

dist_plot <- plot_pca_with_color("Foraging.Mode")

K-Means Clustering

Create Clusters for 2-5 Groups Out of Curiosity

# Create clusters for 2 Groups - To Match Distribution Levels
kmeans_2groups <- data_PCA %>%
  kmeans(centers=2, iter.max = 1000)

# Calculate the within group inertia
# Within Group Inertia
kmeans_2groups$totss

[1] 1005

kmeans_2groups$tot.withinss

[1] 760.1047

kmeans_2groups$withinss

[1] 481.5935 278.5113

kmeans_2groups$betweenss

[1] 244.8953

# Create Clusters For More Groups Out Curiosity
first_kmeans <- data_PCA %>% 
  kmeans(centers = 4, iter.max = 1000)

kmeans_5groups<- data_PCA %>%
  kmeans(centers = 5, iter.max = 1000)

kmeans_4groups<- data_PCA %>%
  kmeans(centers = 4, iter.max = 1000)

kmeans_3groups<- data_PCA %>%
  kmeans(centers=3, iter.max = 1000)

# Extract the Clustering Groups
clusters_5groups <- kmeans_5groups$cluster
clusters_4groups <- kmeans_4groups$cluster
clusters_3groups <- kmeans_3groups$cluster
clusters_2groups <- kmeans_2groups$cluster

Plot a Comparison of the 2 Group K-Means Clustering vs The Coloration Applied by the “Distribution” Variable

The k-means method found clusters which visually seem very similar to those applied by the Distribution category.

Color_label = "Distribution"

# Define a color palette for consistency
custom_palette <- c("red", "blue", "green", "orange")

# Save each of the plots 
p1A <- fviz_pca_biplot(result_pca, axes = c(1, 2), col.ind = as.factor(clusters_2groups), col.var = "black", title = "PC1 & PC2 - Clustering") + ylim(-4, 4)
p1B <- fviz_pca_biplot(result_pca, axes = c(1, 2), col.ind = data_PCA_all_vars[[Color_label]], col.var = "black", title = "PC1 & PC2 - Distribution")+ ylim(-4, 4) 
p2A <- fviz_pca_biplot(result_pca, axes = c(1, 3), col.ind = as.factor(clusters_2groups), col.var = "black", title = "PC1 & PC3 - K-Means Clustering")+ ylim(-4, 4) 
p2B <- fviz_pca_biplot(result_pca, axes = c(1, 3), col.ind = data_PCA_all_vars[[Color_label]], col.var = "black", title = "PC1 & PC3 - Distribution")+ ylim(-4, 4)
p3A <- fviz_pca_biplot(result_pca, axes = c(2, 3), col.ind = as.factor(clusters_2groups), col.var = "black", title = "PC2 & PC3 - K-Means Clustering")+ ylim(-4, 4)
p3B <- fviz_pca_biplot(result_pca, axes = c(2, 3), col.ind = data_PCA_all_vars[[Color_label]], col.var = "black", title = "PC2 & PC3 - Distribution")+ylim(-4, 4)

# Plot each axis combination
all_4_plots<-grid.arrange(p1A, p1B, p2A, p2B, p3A, p3B, nrow = 4, ncol=2)

# Save a ggplot object
ggsave("Clustering-Kmeans.png", plot = all_4_plots, width = 11, height = 10, units = "in", dpi = 500)

Comparison of K-Means Clusterings and Distribution Groupings Using Chi-Squared Test of Independence

# Convert kmeans output to a dataframe
data_with_clusters <- data_PCA %>%
  as.data.frame() %>%                  # Ensure it's a dataframe
  mutate(Cluster = kmeans_2groups$cluster) %>% # Add the cluster assignments
  mutate(Distribution = data_PCA_all_vars$Distribution)
data_with_clusters

# Contingency table
table_kmeans_vs_distribution <- table(kmeans_clusters = data_with_clusters$Cluster, distribution = data_with_clusters$Distribution)

# Print the table
print(table_kmeans_vs_distribution)

               distribution
kmeans_clusters Temperate Tropical
              1        33       93
              2        37       39

# Perform chi-square test
chisq_test <- chisq.test(table_kmeans_vs_distribution)
print(chisq_test)

    Pearson's Chi-squared test with Yates' continuity correction

data:  table_kmeans_vs_distribution
X-squared = 9.6222, df = 1, p-value = 0.001922

Unsupervised Learning

##Split Data ### Specify Variables to Use

# Specify the Variables to Use in the Random Forests Analysis
variables_to_use_RF <- c(
  ### Categorical Variables
  "Mode.of.reproduction",
  "Distribution",
  "Prefered.Habitat.Type",
  "Foraging.Mode",
  
  ### Location
  #"Longitude",
  #"Latitude",
  
  ### Numerical Variables
  #"Average.Female.adult.weight..g.",
  "Mean.F.SVL.adults..mm.",
  #"F.SVL.at.maturity..mm.",
  "Offspring.SVL..mm.",
  "Mean.Clutch.Size",
  "Clutches.per.year",
  #"Clutch.Frequency",
  "RCM"
  
  ### Standard Deviations and Sample Sizes
  #"SD.F.SVL.adults..mm.",
  #"Sample.size.Mean.F.SVL.adults",
  #"SD.Female.adult.weight..g.",
  #"Sample.Size.Female.adult.weight"
)

###Filter Data Variables Filter out the data to only keep variables to be used in the random forests and KNN analyses. Remove Nas.

data_RF <- data %>%
  select(any_of(variables_to_use_RF)) %>% # Select variables
  na.omit() %>%
  droplevels()

Initial Split

Split the data such that 80% is used for the training and 20% is unseen.

split_data <- data_RF %>%
  initial_split(prop = 0.8, strata = "Distribution")

# Train and Test Data
train_data <- analysis(split_data)
test_data <- assessment(split_data)

Random Forests

Choose the Optimum Hyperparameter ‘mtry’.

This chunk enables us to determine the optimum value of the hyperparameter mtry. Here it is found to be at the highest for the mtry values 4, 5 and 6.

# Initialisation
params_CV <- trainControl(method = "repeatedcv", 
                          number = 10, 
                          repeats = 5)

predictor_forest <- caret::train(Distribution ~ ., 
                                 data = data_RF, 
                                 method = "rf", 
                                 metric = "Accuracy", 
                                 trControl = params_CV, 
                                 tuneGrid = expand.grid(mtry = 2:12))

as.data.frame(predictor_forest$results) %>% 
  ggplot(aes(x = mtry, y = Accuracy)) +
  geom_point()

Create the RF

Use 5000 random trees and mtry value 3.

forest <- randomForest(Distribution ~ ., # Formula for prediction
                      data = train_data, # Data for training
                      ntree = 5000, # Number of trees
                      maxnodes = 10, # Number of maximum leaves for each tree
                      mtry = 3, # Number of variables for each tree
                      importance = TRUE) # Computation of importance

Test the Accuracy

Accuracy found to be 0.829

# Fast way of computing accuracy with the pipe operator %>% 
confusion_matrix <-predict(forest, newdata = select(test_data, - Distribution)) %>% 
  table(prediction = ., truth = test_data$Distribution) 

accuracy <- confusion_matrix %>%
  {sum(diag(.)) / sum(.)}

confusion_matrix

           truth
prediction  Temperate Tropical
  Temperate         9        0
  Tropical          5       27

accuracy

[1] 0.8780488

Find the Variable Importance

Found the Mean Clutch Size caused the highest mean decrease in accuracy and mean decrease in the Gini coefficient of impurity. Therefore it is the variable contributing the most to the classification.

png("Variable_importance.png", width = 12, height = 8, units = "in", res = 500)
varImpPlot(
  forest, 
  main = "Variable Importance Plot Random Forests",
  cex = 0.7,
  
)
dev.off()

null device 
          1 

varImpPlot(
  forest, 
  main = "Variable Importance Plot Random Forests",
  cex = 0.7  # Adjust font size
)

KNN

Split the Data Randomly

# We split the data (80% for training)
first_split <- initial_split(data_RF, prop = 0.8)
first_train_set <- analysis(first_split) 
dim(first_train_set); head(first_train_set)

[1] 161   9

first_test_set <- assessment(first_split) 
dim(first_test_set); head(first_test_set)

[1] 41  9

Appendix - CA and CCA

This section is not presented in the main text. We tried to do a CA and CCA but the lack of environmental variables proved to be an obstacle. Absolute latitude was averaged for each habitat in order to test the effect of latitude/climate on the association between families and particular habitats. However no association was found and we were unsure if averaging the value of latitudes for the habitats was appropriate.

Variables to use in the Analysis

variables_to_use_CCA <- c(
  #"Distribution",
  #"Foraging.Mode",
  #"Mode.of.reproduction",
  "Prefered.Habitat.Type",
  "Family",
  "Latitude"
)

Intiial Data Wrangling

# Select all four variables to be used in the CA and CCA, omit NAs
data_CCA <- data %>%
  # Select only the variables to use for the CA and CCA
  select(any_of(variables_to_use_CCA)) %>%
  # Exclude NAs
  na.omit() 

# Make a frequency table containing just 2 categorical variables for the CA
# With each pair of levels on a seperate row
freq_CA_long <- data_CCA %>%
  select(Family, Prefered.Habitat.Type) %>%
  table() %>%
  as_tibble %>%
  rename(Freq = "n",
         Habitat = "Prefered.Habitat.Type")

# Make a tibble with the count of each combination
freq_CA_short <- data_CCA %>%
  select(Prefered.Habitat.Type, Family) %>%
  count(Prefered.Habitat.Type, Family) %>%  
  pivot_wider(
    names_from = Family, 
    values_from = n,  
    values_fill = 0)

# Convert to dataframe and create rownames (for use creating average profile)
freq_CA_short_df <- freq_CA_short %>%
  as.data.frame() %>%  
  column_to_rownames(var = "Prefered.Habitat.Type")  

freq_CA_long

freq_CA_short

freq_CA_short_df

Represent in ggplot

ggplot(freq_CA_long) +
  aes(x = Family, y = Habitat, fill = Freq) + # fill says which column is used 
  # for filling
  geom_raster() + # represents in raster mode
  # and dress to make it pretty
  scale_fill_viridis_c() + # Change color scale
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Calculate the Average Profile

# Create average profile  
average_profile <- colSums(freq_CA_short_df) / sum(freq_CA_short_df)
average_profile

      Dactyloidae  Gymnophthalmidae     Shinisauridae           Teiidae          Agamidae    Chamaeleonidae    Corytophanidae   Diplodactylidae 
      0.054570259       0.065484311       0.001364256       0.169167804       0.036834925       0.002728513       0.002728513       0.015006821 
       Gekkonidae     Hoplocercidae      Leiosauridae       Liolaemidae         Opluridae   Phrynosomatidae  Phyllodactylidae     Polychrotidae 
      0.053206003       0.001364256       0.006821282       0.028649386       0.001364256       0.111869031       0.031377899       0.005457026 
        Scincidae Sphaerodactylidae      Tropiduridae         Varanidae          Anguidae       Anniellidae        Lacertidae        Cordylidae 
      0.163710778       0.025920873       0.077762619       0.004092769       0.010914052       0.004092769       0.038199181       0.006821282 
    Crotaphytidae         Iguanidae    Leiocephalidae      Xenosauridae  Carphodactylidae     Eublepharidae    Helodermatidae       Pygopodidae 
      0.015006821       0.006821282       0.010914052       0.031377899       0.002728513       0.001364256       0.001364256       0.004092769 
      Xantusiidae 
      0.006821282 

# Create the Count  
count_habitat <- rowSums(freq_CA_short_df)
count_habitat

      Aquatic      Arboreal Bromelicolous     Fossorial  Psammophilus    Saxicolous Semi-arboreal   Terrestrial 
            7           133            10            16            10           120             6           431 

# Division Counts
CA_repartition <- (freq_CA_short_df / count_habitat) %>% 
  rbind(Average_profile = average_profile)
CA_repartition

ggplot(freq_CA_long) +
  aes(x = Family, y = Habitat, fill = Freq) + # fill says which column is used 
  # for filling
  geom_raster() + # represents in raster mode
  # and dress to make it pretty
  scale_fill_viridis_c() + # Change color scale
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

CA

# 
unique(data$Family)

 [1] Agamidae          Anguidae          Anniellidae       Carphodactylidae  Chamaeleonidae    Cordylidae        Corytophanidae    Crotaphytidae    
 [9] Dactyloidae       Diplodactylidae   Eublepharidae     Gekkonidae        Gymnophthalmidae  Helodermatidae    Hoplocercidae     Iguanidae        
[17] Lacertidae        Leiocephalidae    Leiosauridae      Liolaemidae       Opluridae         Phrynosomatidae   Phyllodactylidae  Polychrotidae    
[25] Pygopodidae       Scincidae         Shinisauridae     Sphaerodactylidae Teiidae           Tropiduridae      Varanidae         Xantusiidae      
[33] Xenosauridae      <NA>             
33 Levels: Agamidae Anguidae Anniellidae Carphodactylidae Chamaeleonidae Cordylidae Corytophanidae Crotaphytidae Dactyloidae ... Xenosauridae

result_CA <- CA(freq_CA_short_df, graph = FALSE)
fviz_eig(result_CA)

fviz_ca_biplot(result_CA, repel = TRUE)

Canonical Correspondance Analysis

Create Data

mean_latitude_df <- data_CCA %>%
  group_by(Prefered.Habitat.Type) %>%           
  summarise(mean_latitude_abs = mean(abs(Latitude), na.rm = TRUE)) %>% 
  as.data.frame()   

mean_latitude_df

Run CCA

anova(result_CCA,by="margin")

Permutation test for cca under reduced model
Marginal effects of terms
Permutation: free
Number of permutations: 999

Model: cca(formula = freq_CA_short_df ~ mean_latitude_abs, data = mean_latitude_df)
                  Df ChiSquare      F Pr(>F)
mean_latitude_abs  1   0.23944 1.2409  0.761
Residual           6   1.15773